"""
特征提取模块
提取书法作品的结构特征、笔画特征、风格特征等
"""
import cv2
import numpy as np
import pywt
import logging
import traceback
from scipy.fftpack import fft2
from skimage.measure import regionprops
from scipy.spatial import distance
from scipy.stats import variation
from skimage import measure

# 配置日志
logger = logging.getLogger(__name__)

class FeatureExtractor:
    def __init__(self):
        self.features = {}
        self.image = None
        self.binary_image = None
        self.skeleton = None
        self.gray_image = None
        logger.info("初始化特征提取器")
    
    def load_processed_images(self, processed_images):
        """加载预处理后的图像"""
        try:
            logger.info("加载预处理图像")
            self.image = processed_images.get('original')
            self.binary_image = processed_images.get('binary')
            self.skeleton = processed_images.get('skeleton')
            self.gray_image = processed_images.get('grayscale')
            
            if self.binary_image is None:
                logger.error("二值图像不存在")
                raise ValueError("二值图像不存在，请先进行图像预处理")
            
            # 记录加载的图像信息
            if self.image is not None:
                logger.debug(f"原始图像形状: {self.image.shape}, 类型: {self.image.dtype}")
            if self.binary_image is not None:
                logger.debug(f"二值图像形状: {self.binary_image.shape}, 类型: {self.binary_image.dtype}")
            if self.skeleton is not None:
                logger.debug(f"骨架图像形状: {self.skeleton.shape}, 类型: {self.skeleton.dtype}")
            if self.gray_image is not None:
                logger.debug(f"灰度图像形状: {self.gray_image.shape}, 类型: {self.gray_image.dtype}")
                
            logger.info("图像加载完成")
            return self
        except Exception as e:
            logger.error(f"加载图像时出错: {str(e)}")
            logger.error(traceback.format_exc())
            raise
        
    def extract_all_features(self):
        """提取所有特征"""
        try:
            logger.info("开始提取所有特征")
            
            # 空间结构特征
            logger.info("提取空间结构特征")
            self.extract_spatial_features()
            
            # 笔画动态特征
            logger.info("提取笔画动态特征")
            self.extract_stroke_features()
            
            # 风格特征
            logger.info("提取风格特征")
            self.extract_style_features()
            
            # 整篇作品章法特征 
            logger.info("提取整篇作品章法特征")
            self.extract_composition_features()
            
            # 记录提取的特征信息
            logger.info(f"特征提取完成，共{len(self.features)}个特征")
            for key, value in self.features.items():
                try:
                    if isinstance(value, list):
                        logger.debug(f"特征 {key}: 类型={type(value)}, 长度={len(value)}, "
                                    f"元素类型={type(value[0]) if value else 'empty'}")
                    elif hasattr(value, 'shape'):
                        logger.debug(f"特征 {key}: 类型={type(value)}, 形状={value.shape}, 数据类型={value.dtype}")
                    else:
                        logger.debug(f"特征 {key}: 类型={type(value)}, 值={value}")
                except Exception as debug_error:
                    logger.debug(f"无法记录特征 {key} 的信息: {str(debug_error)}")
            
            return self.features
        except Exception as e:
            logger.error(f"提取特征时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return self.features  # 返回已提取的特征
    
    def extract_spatial_features(self):
        """提取空间结构特征"""
        # 1. 九宫格区域密度分布
        self.features['grid_density'] = self._calculate_grid_density()
        
        # 2. 字体重心偏移量
        self.features['center_offset'] = self._calculate_center_offset()
        
        # 3. 部件密度分布
        self.features['component_distribution'] = self._calculate_component_distribution()
        
        return self.features
    
    def extract_stroke_features(self):
        """提取笔画动态特征"""
        # 1. 线条曲率变化
        self.features['curvature'] = self._calculate_curvature()
        
        # 2. 角度分布
        self.features['angle_distribution'] = self._calculate_angle_distribution()
        
        # 3. 线条宽度变化
        self.features['line_width_variance'] = self._calculate_line_width_variance()
        
        return self.features
    
    def extract_style_features(self):
        """提取风格特征"""
        # 1. 纹理特征（小波域能量分布）
        self.features['wavelet_energy'] = self._calculate_wavelet_energy()
        
        # 2. 形状特征（傅里叶描述子）
        self.features['fourier_descriptors'] = self._calculate_fourier_descriptors()
        
        # 3. 整体结构特征
        self.features['structural_complexity'] = self._calculate_structural_complexity()
        
        return self.features
    
    def extract_composition_features(self):
        """提取整篇作品章法特征"""
        # 1. 平衡度 - 作品整体视觉重心的平衡性
        self.features['visual_balance'] = self._calculate_visual_balance()
        
        # 2. 行列间距一致性 - 评估多字作品中行距和列距的一致性
        self.features['spacing_consistency'] = self._calculate_spacing_consistency()
        
        # 3. 整体空间利用率 - 作品充分利用空间的程度
        self.features['space_utilization'] = self._calculate_space_utilization()
        
        # 4. 边缘留白比例 - 作品四周留白的均衡性
        self.features['margin_ratio'] = self._calculate_margin_ratio()
        
        return self.features
    
    def _calculate_grid_density(self, grid_size=3):
        """计算九宫格区域密度分布"""
        if self.binary_image is None:
            return None
        
        h, w = self.binary_image.shape
        grid_h, grid_w = h // grid_size, w // grid_size
        
        densities = []
        
        # 黑色为笔画（值为0）
        for i in range(grid_size):
            for j in range(grid_size):
                cell = self.binary_image[i*grid_h:(i+1)*grid_h, j*grid_w:(j+1)*grid_w]
                # 计算笔画像素的比例
                density = float(np.sum(cell == 0) / (grid_h * grid_w))
                densities.append(density)
        
        return densities
    
    def _calculate_center_offset(self):
        """计算书法中心偏移特征"""
        try:
            logger.debug("开始计算中心偏移")
            
            if self.binary_image is None:
                logger.warning("二值图像不存在，无法计算中心偏移")
                return [0.0, 0.0]
            
            # 找到所有笔画像素点
            logger.debug("查找笔画像素点")
            stroke_pixels = np.where(self.binary_image > 0)
            logger.debug(f"找到{len(stroke_pixels[0])}个笔画像素点")
            
            if len(stroke_pixels[0]) == 0:
                logger.warning("没有找到笔画像素点")
                return [0.0, 0.0]
            
            # 计算图像的几何中心
            logger.debug("计算图像几何中心")
            img_height, img_width = self.binary_image.shape
            img_center_x = img_width / 2.0
            img_center_y = img_height / 2.0
            logger.debug(f"图像几何中心: ({img_center_x}, {img_center_y})")
            
            # 计算笔画的质心
            logger.debug("计算笔画质心")
            stroke_center_y = float(np.mean(stroke_pixels[0]))
            stroke_center_x = float(np.mean(stroke_pixels[1]))
            logger.debug(f"笔画质心: ({stroke_center_x}, {stroke_center_y})")
            
            # 计算相对偏移（归一化到[-1, 1]范围）
            logger.debug("计算相对偏移")
            x_offset = float((stroke_center_x - img_center_x) / (img_width / 2.0))
            y_offset = float((stroke_center_y - img_center_y) / (img_height / 2.0))
            logger.debug(f"X偏移: {x_offset}, Y偏移: {y_offset}")
            
            # 返回列表而不是元组
            return [float(x_offset), float(y_offset)]
        except Exception as e:
            logger.error(f"计算中心偏移时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return [0.0, 0.0]
    
    def _calculate_component_distribution(self):
        """计算组件分布特征"""
        try:
            logger.debug("开始计算组件分布特征")
            
            if self.binary_image is None:
                logger.warning("二值图像不存在，无法计算组件分布特征")
                return [0.0, 0.0, 0.0, 0.0]
            
            # 查找连通区域
            logger.debug("查找连通区域")
            num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(self.binary_image, connectivity=8)
            logger.debug(f"找到{num_labels}个连通区域")
            
            # 忽略背景（通常是第一个）
            stats = stats[1:]
            centroids = centroids[1:]
            
            # 如果没有找到有效的连通区域或只找到一个（可能是背景）
            num_features = len(stats)
            logger.debug(f"有效连通区域数量: {num_features}")
            
            if num_features <= 1:
                logger.warning("连通区域数量不足，无法计算有意义的分布特征")
                return [0.0, 0.0, 0.0, 0.0]
            
            # 计算组件中心点的分布特征
            logger.debug("计算中心点坐标")
            x_coords = [centroid[0] for centroid in centroids]
            y_coords = [centroid[1] for centroid in centroids]
            
            # 计算标准差、方差系数等特征
            logger.debug("计算统计特征")
            x_std = float(np.std(x_coords))
            y_std = float(np.std(y_coords))
            
            # 计算变异系数
            x_mean = float(np.mean(x_coords))
            y_mean = float(np.mean(y_coords))
            
            x_cv = float(x_std / x_mean) if x_mean != 0 else 0.0
            y_cv = float(y_std / y_mean) if y_mean != 0 else 0.0
            
            logger.debug(f"X坐标标准差: {x_std}, Y坐标标准差: {y_std}")
            logger.debug(f"X坐标变异系数: {x_cv}, Y坐标变异系数: {y_cv}")
            
            # 返回四个基本特征
            return [float(x_std), float(y_std), float(x_cv), float(y_cv)]
        except Exception as e:
            logger.error(f"计算组件分布特征时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return [0.0, 0.0, 0.0, 0.0]
    
    def _calculate_curvature(self):
        """计算书法曲率特征"""
        try:
            logger.debug("开始计算曲率特征")
            
            if self.skeleton is None:
                logger.warning("骨架图像不存在，无法计算曲率特征")
                return [0.0, 0.0, 0.0]
            
            # 查找骨架中的点
            logger.debug("查找骨架点")
            y_coords, x_coords = np.where(self.skeleton == 0)
            points = np.column_stack((x_coords, y_coords))
            logger.debug(f"找到{len(points)}个骨架点")
            
            if len(points) < 3:
                logger.warning("骨架点数量不足以计算曲率")
                return [0.0, 0.0, 0.0]
            
            # 计算曲率特征
            try:
                logger.debug("计算每个点的曲率")
                curvatures = []
                
                # 遍历每个骨架点，计算局部曲率
                for i in range(1, len(points) - 1):
                    # 获取相邻三个点
                    p_prev = points[i-1]
                    p_curr = points[i]
                    p_next = points[i+1]
                    
                    # 计算向量
                    vec1 = p_curr - p_prev
                    vec2 = p_next - p_curr
                    
                    # 计算向量长度
                    len1 = np.linalg.norm(vec1)
                    len2 = np.linalg.norm(vec2)
                    
                    # 避免除零错误
                    if len1 > 0 and len2 > 0:
                        # 归一化向量
                        vec1_norm = vec1 / len1
                        vec2_norm = vec2 / len2
                        
                        # 计算点积，得到夹角的余弦值
                        dot_product = np.dot(vec1_norm, vec2_norm)
                        
                        # 确保dot_product在[-1, 1]范围内，避免数值误差
                        dot_product = np.clip(dot_product, -1.0, 1.0)
                        
                        # 计算角度
                        angle = np.arccos(dot_product)
                        
                        # 计算曲率(夹角/弧长)
                        curvature = angle / ((len1 + len2) / 2)
                        curvatures.append(curvature)
                
                logger.debug(f"计算得到{len(curvatures)}个曲率值")
                
                if not curvatures:
                    logger.warning("没有计算出曲率值")
                    return [0.0, 0.0, 0.0]
                
                # 计算统计特征
                curvature_mean = float(np.mean(curvatures))
                curvature_std = float(np.std(curvatures))
                curvature_max = float(np.max(curvatures))
                
                logger.debug(f"曲率均值: {curvature_mean}, 标准差: {curvature_std}, 最大值: {curvature_max}")
                return [curvature_mean, curvature_std, curvature_max]
            
            except Exception as inner_e:
                logger.error(f"计算曲率过程中出错: {str(inner_e)}")
                logger.error(traceback.format_exc())
                return [0.0, 0.0, 0.0]
                
        except Exception as e:
            logger.error(f"计算曲率特征时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return [0.0, 0.0, 0.0]
    
    def _calculate_angle_distribution(self, num_bins=8):
        """计算角度分布直方图"""
        try:
            logger.debug("开始计算角度分布")
            
            if self.skeleton is None:
                logger.warning("骨架图像不存在，无法计算角度分布")
                return None
            
            # 使用Sobel算子计算梯度
            grad_x = cv2.Sobel(self.skeleton, cv2.CV_64F, 1, 0, ksize=3)
            grad_y = cv2.Sobel(self.skeleton, cv2.CV_64F, 0, 1, ksize=3)
            logger.debug(f"梯度X形状: {grad_x.shape}, 梯度Y形状: {grad_y.shape}")
            
            # 计算角度和幅值
            mag, angle = cv2.cartToPolar(grad_x, grad_y)
            logger.debug(f"角度形状: {angle.shape}, 幅值形状: {mag.shape}")
            
            # 只考虑骨架上的像素点（值为0的点）
            valid_indices = np.where(self.skeleton == 0)
            valid_angles = angle[valid_indices]
            valid_mags = mag[valid_indices]
            logger.debug(f"有效角度数量: {len(valid_angles)}")
            
            if len(valid_angles) == 0:
                logger.warning("没有有效的角度值")
                return [0.0] * num_bins
            
            # 将角度范围（0-2π）划分为num_bins个等分，计算直方图
            logger.debug("计算角度直方图")
            hist, bin_edges = np.histogram(valid_angles, bins=num_bins, range=(0, 2*np.pi), weights=valid_mags)
            logger.debug(f"直方图形状: {hist.shape}, 类型: {hist.dtype}, 分桶边界: {bin_edges}")
            
            # 归一化直方图
            hist_sum = np.sum(hist)
            if hist_sum > 0:
                logger.debug(f"直方图和: {hist_sum}")
                hist = hist / hist_sum
            else:
                logger.warning("直方图和为零")
            
            # 确保返回的是Python基本数据类型列表，而不是NumPy数组
            result = [float(val) for val in hist]
            logger.debug(f"角度分布结果: {result}")
            return result
        except Exception as e:
            logger.error(f"计算角度分布时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return [0.0] * num_bins
    
    def _calculate_line_width_variance(self):
        """计算笔画宽度的变化程度"""
        try:
            logger.debug("开始计算笔画宽度变化")
            
            if self.binary_image is None:
                logger.warning("二值图像不存在，无法计算笔画宽度变化")
                return 0.0
                
            # 使用距离变换计算每个像素到背景的距离
            logger.debug("执行距离变换")
            dist_transform = cv2.distanceTransform(self.binary_image, cv2.DIST_L2, 5)
            logger.debug(f"距离变换结果形状: {dist_transform.shape}, 类型: {dist_transform.dtype}")
            
            # 找到所有笔画像素点（前景）
            stroke_pixels = dist_transform[self.binary_image > 0]
            logger.debug(f"笔画像素点数量: {len(stroke_pixels)}")
            
            if len(stroke_pixels) == 0:
                logger.warning("没有找到笔画像素点")
                return 0.0
                
            # 计算变异系数 (CV)
            logger.debug("计算变异系数")
            cv_value = float(variation(stroke_pixels))
            
            # 处理变异系数可能为NaN的情况
            if np.isnan(cv_value):
                logger.warning("变异系数为NaN")
                return 0.0
                
            logger.debug(f"笔画宽度变异系数: {cv_value}")
            return float(cv_value)
        except Exception as e:
            logger.error(f"计算笔画宽度变化时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return 0.0
    
    def _calculate_wavelet_energy(self, wavelet='db1', level=3):
        """计算小波域纹理能量分布"""
        if self.gray_image is None:
            if hasattr(self, 'image') and self.image is not None:
                if len(self.image.shape) == 3:
                    gray = cv2.cvtColor(self.image, cv2.COLOR_BGR2GRAY)
                else:
                    gray = self.image
            else:
                return [0.0] * (1 + 3 * level)  # 返回全零向量
        else:
            gray = self.gray_image
        
        # 调整图像大小为2的幂次，便于小波分解
        h, w = gray.shape
        new_h, new_w = 2**int(np.log2(h)), 2**int(np.log2(w))
        gray_resized = cv2.resize(gray, (new_w, new_h))
        
        # 二维小波分解
        coeffs = pywt.wavedec2(gray_resized, wavelet, level=level)
        
        # 计算每个子带的能量
        energy_features = []
        
        # 添加低频子带能量
        cA = coeffs[0]
        energy_A = float(np.sum(cA**2) / cA.size)
        energy_features.append(energy_A)
        
        # 添加高频子带能量
        for i in range(1, len(coeffs)):
            (cH, cV, cD) = coeffs[i]
            energy_H = float(np.sum(cH**2) / cH.size)
            energy_V = float(np.sum(cV**2) / cV.size)
            energy_D = float(np.sum(cD**2) / cD.size)
            energy_features.extend([energy_H, energy_V, energy_D])
        
        # 归一化能量特征
        total_energy = np.sum(energy_features)
        if total_energy > 0:
            energy_features = [float(e / total_energy) for e in energy_features]
        
        return energy_features
    
    def _calculate_fourier_descriptors(self, num_descriptors=20):
        """计算轮廓的傅里叶描述子"""
        try:
            logger.debug("开始计算傅里叶描述子")
            
            if self.binary_image is None:
                logger.warning("二值图像不存在，无法计算傅里叶描述子")
                return [0.0] * num_descriptors
            
            # 提取轮廓
            logger.debug("提取轮廓")
            contours, _ = cv2.findContours(255 - self.binary_image, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
            logger.debug(f"找到{len(contours)}个轮廓")
            
            # 如果没有找到轮廓，返回零向量
            if len(contours) == 0:
                logger.warning("没有找到轮廓")
                return [0.0] * num_descriptors
            
            # 找到最大轮廓
            max_contour = max(contours, key=cv2.contourArea)
            logger.debug(f"最大轮廓点数: {len(max_contour)}")
            
            # 如果轮廓点数太少，返回零向量
            if len(max_contour) < num_descriptors:
                logger.warning(f"轮廓点数({len(max_contour)})小于所需描述子数({num_descriptors})")
                return [0.0] * num_descriptors
            
            # 将轮廓转换为复数表示形式
            logger.debug("将轮廓转换为复数表示")
            contour_complex = np.empty(len(max_contour), dtype=complex)
            for i in range(len(max_contour)):
                contour_complex[i] = complex(max_contour[i][0][0], max_contour[i][0][1])
            
            # 计算傅里叶变换
            logger.debug("计算傅里叶变换")
            fourier_result = np.fft.fft(contour_complex)
            logger.debug(f"傅里叶结果形状: {fourier_result.shape}, 类型: {fourier_result.dtype}")
            
            # 计算幅度谱
            fourier_amplitudes = np.absolute(fourier_result)
            logger.debug(f"幅度谱形状: {fourier_amplitudes.shape}")
            
            # 归一化
            logger.debug("归一化幅度谱")
            fourier_amplitudes[0] = 0  # 忽略直流分量
            norm_factor = np.sum(fourier_amplitudes)
            
            # 避免除零错误
            if float(norm_factor) > 0:
                logger.debug(f"归一化因子: {norm_factor}")
                normalized_amplitudes = fourier_amplitudes / float(norm_factor)
            else:
                logger.warning("归一化因子为零")
                return [0.0] * num_descriptors
            
            # 取前num_descriptors个描述子
            descriptors = list(normalized_amplitudes[:num_descriptors])
            
            # 确保返回的是Python基本数据类型，而不是NumPy数组元素
            result = [float(val) for val in descriptors]
            logger.debug(f"傅里叶描述子: {result[:5]}... (共{len(result)}个)")
            return result
        except Exception as e:
            logger.error(f"计算傅里叶描述子时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return [0.0] * num_descriptors
    
    def _calculate_structural_complexity(self):
        """计算结构复杂度"""
        try:
            logger.debug("开始计算结构复杂度")
            
            if self.binary_image is None:
                logger.warning("二值图像不存在，无法计算结构复杂度")
                return 0.0
                
            # 计算轮廓
            logger.debug("查找轮廓")
            contours, _ = cv2.findContours(self.binary_image, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
            logger.debug(f"找到{len(contours)}个轮廓")
            
            if not contours:
                logger.warning("没有找到轮廓")
                return 0.0
                
            # 计算总周长和总面积
            total_perimeter = 0
            total_area = 0
            
            for cont in contours:
                perimeter = cv2.arcLength(cont, True)
                area = cv2.contourArea(cont)
                logger.debug(f"轮廓周长: {perimeter}, 面积: {area}")
                
                if area > 0:  # 避免除零错误
                    total_perimeter += perimeter
                    total_area += area
            
            if total_area == 0:
                logger.warning("总面积为零")
                return 0.0
                
            # 计算复杂度 (周长^2 / 面积)
            complexity = (total_perimeter ** 2) / total_area
            logger.debug(f"结构复杂度: {complexity}")
            
            # 归一化复杂度 (使用sigmoid函数限制在[0,1]范围)
            normalized_complexity = 1 / (1 + np.exp(-complexity / 1000 + 5))
            logger.debug(f"归一化后的结构复杂度: {normalized_complexity}")
            
            return float(normalized_complexity)
        except Exception as e:
            logger.error(f"计算结构复杂度时出错: {str(e)}")
            logger.error(traceback.format_exc())
            return 0.0
    
    def _calculate_visual_balance(self):
        """计算作品视觉平衡度
        
        视觉平衡度衡量作品的重心分布是否均衡，理想的作品应在视觉上保持平衡
        """
        if self.binary_image is None:
            return 0.5  # 默认中等平衡度
        
        # 计算图像重心
        h, w = self.binary_image.shape
        y_indices, x_indices = np.where(self.binary_image > 0)
        
        if len(y_indices) == 0:
            return 0.5
        
        center_y = np.mean(y_indices)
        center_x = np.mean(x_indices)
        
        # 计算图像中心点
        img_center_y = h / 2
        img_center_x = w / 2
        
        # 计算重心偏移比例
        offset_y = abs(center_y - img_center_y) / (h / 2)
        offset_x = abs(center_x - img_center_x) / (w / 2)
        
        # 归一化为0-1之间的平衡度指标，1表示完美平衡
        balance = 1 - (offset_x + offset_y) / 2
        
        logger.debug(f"计算视觉平衡度: 重心=({center_x:.2f}, {center_y:.2f}), 偏移=({offset_x:.4f}, {offset_y:.4f}), 平衡度={balance:.4f}")
        return balance
    
    def _calculate_spacing_consistency(self):
        """计算行列间距一致性
        
        对于多字作品，评估行距和列距的一致性，检测笔画或部件之间的间距是否均匀。
        高一致性意味着整齐有序的排列，有助于形成良好的章法。
        """
        if self.binary_image is None:
            return [0.5, 0.5]  # 默认中等一致性
        
        try:
            # 使用连通区域分析查找可能的字符区域
            labeled_img = measure.label(self.binary_image)
            regions = measure.regionprops(labeled_img)
            
            # 如果区域太少，可能是单字作品
            if len(regions) < 3:
                logger.debug("连通区域太少，可能不是多字作品，返回默认间距一致性")
                return [0.7, 0.7]
            
            # 提取所有区域的中心点
            centers = np.array([region.centroid for region in regions])
            
            # 计算水平和垂直方向上的距离
            if len(centers) > 1:
                # 排序中心点按y坐标（行）
                sorted_by_y = centers[np.argsort(centers[:, 0])]
                
                # 计算相邻行之间的距离
                row_dists = np.diff(sorted_by_y[:, 0])
                row_consistency = 1.0
                if len(row_dists) > 0:
                    # 计算行距变异系数 (标准差/均值)，越小越一致
                    cv_row = np.std(row_dists) / np.mean(row_dists) if np.mean(row_dists) > 0 else 0
                    row_consistency = max(0, 1 - min(cv_row, 1))
                
                # 排序中心点按x坐标（列）
                sorted_by_x = centers[np.argsort(centers[:, 1])]
                
                # 计算相邻列之间的距离
                col_dists = np.diff(sorted_by_x[:, 1])
                col_consistency = 1.0
                if len(col_dists) > 0:
                    # 计算列距变异系数，越小越一致
                    cv_col = np.std(col_dists) / np.mean(col_dists) if np.mean(col_dists) > 0 else 0
                    col_consistency = max(0, 1 - min(cv_col, 1))
                
                logger.debug(f"行距一致性: {row_consistency:.4f}, 列距一致性: {col_consistency:.4f}")
                return [row_consistency, col_consistency]
            else:
                logger.debug("无法检测到多个中心点，返回默认间距一致性")
                return [0.7, 0.7]
        except Exception as e:
            logger.warning(f"计算间距一致性时出错: {str(e)}")
            return [0.5, 0.5]
    
    def _calculate_space_utilization(self):
        """计算整体空间利用率
        
        评估书法作品对整体空间的利用程度，良好的章法应该有合理的空间利用率，
        既不过于拥挤也不过于稀疏。
        """
        if self.binary_image is None:
            return 0.5
        
        # 计算墨迹像素占整体像素的比例
        h, w = self.binary_image.shape
        total_pixels = h * w
        ink_pixels = np.sum(self.binary_image > 0)
        
        # 直接使用墨迹比例作为特征值可能不合适，因为理想的利用率存在最优区间
        ink_ratio = ink_pixels / total_pixels
        
        # 转换为0-1的评分，假设最佳利用率在15%-25%之间
        # 过低或过高都会降低评分
        if 0.15 <= ink_ratio <= 0.25:
            score = 1.0
        elif ink_ratio < 0.15:
            # 墨迹太少，空间利用不足
            score = ink_ratio / 0.15
        else:  # ink_ratio > 0.25
            # 墨迹太多，过于拥挤
            score = max(0, 1 - (ink_ratio - 0.25) / 0.25)
            
        logger.debug(f"空间利用率: 墨迹比例={ink_ratio:.4f}, 评分={score:.4f}")
        return score
    
    def _calculate_margin_ratio(self):
        """计算边缘留白比例
        
        分析作品四周留白的分布情况，理想的章法应有适当且均衡的留白。
        返回四个边缘(上、右、下、左)的留白比例。
        """
        if self.binary_image is None:
            return [0.2, 0.2, 0.2, 0.2]  # 默认均匀留白
        
        h, w = self.binary_image.shape
        
        # 找到墨迹的边界框
        y_indices, x_indices = np.where(self.binary_image > 0)
        
        if len(y_indices) == 0:
            return [0.2, 0.2, 0.2, 0.2]
        
        top = np.min(y_indices)
        bottom = np.max(y_indices)
        left = np.min(x_indices)
        right = np.max(x_indices)
        
        # 计算四个方向的留白比例
        top_margin = top / h
        bottom_margin = (h - bottom - 1) / h
        left_margin = left / w
        right_margin = (w - right - 1) / w
        
        margins = [top_margin, right_margin, bottom_margin, left_margin]
        
        logger.debug(f"留白比例: 上={top_margin:.4f}, 右={right_margin:.4f}, "
                    f"下={bottom_margin:.4f}, 左={left_margin:.4f}")
        
        return margins 